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The development and verification of the Charring Ablating Thermal Protection Implicit System Solver 
(CATPISS) is presented. This work concentrates on the derivation and verification of the stationary grid terms 
in the equations that govern three-dimensional heat and mass transfer for charring thermal protection systems 
including pyrolysis gas flow through the porous char layer. The governing equations are discretized according 
to the Galerkin finite element method (FEM) with first and second order fully implicit time integrators. The 
governing equations are fully coupled and are solved in parallel via Newton’s method, while the linear system 
is solved via the Generalized Minimum Residual method (GMRES). Verification results from exact solutions 
and Method of Manufactured Solutions (MMS) are presented to show spatial and temporal orders of accuracy 
as well as nonlinear convergence rates. 


I. Introduction 


Cite Amar 1 


II. Mathematical Model 


A. Governing Equations 


The equations that govern the solid/gas system of the porous charring ablator include energy and mass conservation 
equations for the solid as well as the Navier-Stokes equations as applied to all of the gaseous species considered. In 
the general case, it is possible that the pyrolysis gases react with the remaining solid, or deposit residue (coke) on 
the solid, but these phenomena are neglected. Under the assumptions that the pyrolysis gas is in thermochemical 
equilibrium and the solid and gas are in thermal equilibrium, then the solid and gas energy equations for a stationary 
grid reduce to a mixture energy equation given by 

= V • (kVr) - V • {tPgho'Vg) + Q (1) 


where p, e 0 , </>, h 0 , v, and Q denote density, total energy, porosity, total enthalpy, velocity, and volumetric energy 
source respectively, and the subscript g denotes a quantity with respect to the pyrolysis gases. The velocity of the 
pyrolysis gases is governed by a porous flow law such as Darcy’s law. Since ablators in general can be anisotropic 
materials, the thermal conductivity, k, is a second order tensor. The solid mass conservation equation is simply 


dp s 

dt 


= m s 


( 2 ) 


If it is assumed that all solid decomposition results in pyrolysis gas generation, the gases are free to flow through the 
porous medium, and the gases occupy all of the pore space, then the gas mass conservation equation is given by 


dQp g ) = 

dt 




( 3 ) 
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Before manipulating the governing equations to yield a form suitable for implementation in a finite-element frame- 
work, it is first necessary to discuss the material model that will characterize the two-phase system. 

B. Material Model 

In order to sufficiently explain the governing equations and boundary conditions, it is important to understand the 
material model used to characterize the state of the solid/gas mixture. It is assumed that all the pores are intercon- 
nected, and therefore pyrolysis gases occupy all of the pore space and are free to flow through it. Consequently, the 
density of the solid/gas mixture is described by 


P=<t>Pg + Pa ( 4 ) 

where the solid density is a bulk density, the gas density is a density with respect to the space the gas occupies (pore 
space), and the porosity is equal to the gas volume fraction. In terms of units, Eq. 4 can be expressed as 


P 4> Pg Ps 

[total mass] [pore vol] [gas mass] [solid mass] 
[total vol] [total vol] [pore vol] [total vol] 


( 5 ) 


It is assumed that the thermodynamic state of the pyrolysis gases can be described as a mixture of perfect gases, and 
that the solid and gas phases are in thermal equilibrium resulting in 


T g = T s = T ( 6 ) 

P = f(p g ,T) (7) 

where the equation of state and all relationships describing the thermochemical state of the gas are dictated by an 
external equilibrium chemistry module such as Cantera, CEA, MUTATION, etc. During the code development effort, 
Cantera has been used to provide the thermochemical state of the gas. 

The solid material model adopted in this study is similar to the model developed by Moyer and Rindal (need 
citation) but has been expanded to include an arbitrary number of components, nc. The solid bulk density is given by 


Ps = ^2 TlPi ( 8 ) 

i=l 


where k, is the volume fraction of the i th component in the virgin composite and is therefore constant. The units 
associated with the solid bulk density model in Eq. 8 are 


Pa 

[solid mass] 
[total vol] 




?:= t 


[initial vol of i th comp.] [mass of i th comp.] 


[total vol] 


[initial vol of i th comp.] 


( 9 ) 


It is assumed that the solid does not change volume due to thermal expansion, and therefore the total volume is 
constant. It is important to note that the solid description in Eq. 8 is only a modeling assumption, and the solid is not 
truly comprised of nc components, species, or distinguishable materials. This modeling assumption comes as a result 
of decomposition data obtained from thermogravimetric analyses (TGA). It has been observed that phenolic resins 
undergo a two-stage decomposition process that can be appropriately captured by a two resin component model (need 
citation). 

It is assumed that all decomposed solid mass results in gas mass generation, and the general model of the decom- 
position process is described by 

virgin plastic => char + gas 


This is a generalized description of the initial and final states of the system between which there are transitional states. 
The reaction is irreversible, and the pyrolysis gases are assumed to be in thermochemical equilibrium and do not react 
with the remaining solid in the pore space. 

Taking the temporal derivative of Eq. 8 gives the solid decomposition rate in terms of component decomposition 
rates. 


dp s __ \ ' p dfh 
8t ~ ^ 1 dt 

i—1 


( 10 ) 
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It is assumed that the decomposition of each component can be described by an Arrhenius relationship of the form 


dpt 

at 


kiPvi 



4>i 

e ~ E i/RT for ■ = B and c 


(ID 


which applies at a constant spatial location (as apposed to a given node which can move during the solution process). 

Since most thermophysical properties of the solid are only known for the virgin plastic and fully charred states, 
the intermediate solid is modeled as some interpolated state between virgin and char. This interpolated state is 
characterized by the extent of reaction (/?), or degree of char, given by 


P Pv Ps 

Pv Pc 


( 12 ) 


where the virgin and char bulk densities are known constants. It is evident that as the solid decomposes from virgin 
to char, the extent of reaction ranges from 0 to 1. The definition in Eq. 12 can be rearranged to more clearly describe 
the interpolated state. 

p s = (1 - /3) p v + (3p c (13) 

Although the virgin and char materials are not distinguishable entities within the intermediate solid, Eq. 13 reveals 
that the degree of char represents an effective char volume fraction within the solid (not in the solid/gas mixture). In 
a similar light, CMA defines an effective virgin mass fraction given by 


Pv 


Vv = 

Pv Pc 

which can be related to the extent of reaction through 

Pv 


1-* 


Similarly the char mass fraction is given by 


!/„ = -( 1-/3) 

Ps 


1 rc r, 

y c =i--y v = —p 

Ps 


These effective parameters are used to determine several solid and mixture properties. 


(14) 


(15) 

(16) 


C. Porous Flow Laws 

Applying the Navier-Stokes momentum equations to flow through the char layer would require detailed knowledge 
of the pore structure, and that information is typically not known. Consequently, a porous flow law can be used as 
a simplified momentum equation. Porous flow laws typically require extra knowledge about the material beyond 
thermophysical properties. These properties include the porosity, <f>, and permeability, k of the solid, as well as the 
viscosity, //, of the gas flowing through the porous medium. The porosity and permeability can be determined through 
material testing and is provided to CATPISS as a function of extent of reaction, /3, which is given by Eq. 12. 


1. Darcy’s Law 


In 1856, Darcy (cite Darcy) published results from a series of experiments in which he determined how the vol- 
umetric flow rate, Q, of a laminar flowing fluid relates to the local pressure gradient within a fully saturated porous 
medium. _ 

Q = -A-VP (17) 

P 

where k, is the anisotropic permeability tensor. The superficial or filtration velocity is the volumetric flow rate averaged 
over the cross-sectional area of the medium and is given by 


-i=r7 VP (18) 

The average or seepage velocity of the fluid is the volumetric flow rate averaged over the cross-sectional area through 
which the fluid can flow (porous area) and is given by 


v 


9 


Q_ 

4>A 


-fvp 

< pp 


( 19 ) 
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which assumes that the surface porosity is equal to the volumetric porosity. Darcy’s law is valid for steady laminar 
flows with ’’sufficiently” low Reynolds numbers. While is is conceded that the pyrolysis gas flow within the porous 
ablator can not be adequately described by Darcy’s law, it is being used for code development purposes due to its 
simplicity. 

The seepage velocity can be use to determine the gas mass flux at any point within the medium according to 

trig = (<? bp g ) v ff = - (<f)p g ) A VP (20) 

Moving all terms to the LHS and simplifying gives 

tfiflVg + P = 0 (21) 

For the k th dimension we have 

( v g)k + • VP = 0 (22) 

where Kp- denotes the k th row vector of the permeability tensor. 


D. Property Models 

1. Internal Energy and Enthalpy 

The total internal energy and of the system can be described by 

pe 0 — (1 P ) Pv^-v T fipc&c T (/>P g e 0g 

where 

e v/c = K /c = h° v/c + [ C PWc (TO dT ’ 

JT 0 

and the total energy of the gas is 

e °s = e °g + J T Cv g (T 1 ) dT 1 + - (v g ■ Vg) 

Similarly, the total enthalpy of the gas is 

K = h ° 9 + £ °p a ( r ) dr + \ (v fl ■ v fl ) 


(23) 

(24) 


(25) 


(26) 


2. Thermal Conductivity 


The thermal conductivity is in general an anisotropic function of temperature. The thermal conductivity tensor is 
given by 


k(T) 


fcn(T) k \2 (T) ki 3 (T) 
*21 (T) fc 22 (T) fc 23 (T) 

*31 (T) *32 (T) *33 (T) 


(27) 


where each component of the tensor can be an independent function of temperature. For materials with anisotropies 
not aligned with coordinate axes, an additional transformation matrix must be introduced to appropriately apply the 
model. The thermal conductivity of the mixture is assumed to be a mass weighted average of virgin char and gas. 


* = 


(1 - P) p v 
P 



(28) 


3. Permeability 


The permeability is in general an anisotropic function of the extent of reaction defined by Eq. 12. The permeability 
tensor is given by 


k(P) 


Ull (P) «12 (/?) «13 (P) 

«21 (P) n 2 2 (P) K23 (P) 
«31 (P) n 32 (P) k 33 (P) 


( 29 ) 
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where each component of the tensor can be an independent function of extent of reaction. For materials with 
anisotropies not aligned with coordinate axes, an additional transformation matrix must be introduced to appropri- 
ately apply the model. 

4. Porosity 

The porosity is in general a function of extend of reaction defined by Eq. 12 and represents the gas mass fraction 
since all of the pore space is assumed to be interconnected. 

<!> = HP) (30) 


5. Emissivity and Absorptivity 

The emissivity and absorptivity of the mixture are assumed to be equal to the emissivity and absorptivity of the 
solid. They are modeled similarly to the thermal conductivity while only taking the solid portion into account 


and 


(1 - (3) Pv , (3p c 
e = e v H e c 

Ps Ps 


(31) 


(1 - (3) p v f3p c 

ol — ql v a c 

Ps Ps 


(32) 


E. Galerkin Weak Statement 

Given the material model, the governing equations can be manipulated to give 


Energy: 


d(j> de, 
' 9 dTdp g 


2 y-fyy w^ g 

pCy Pq nrr, O . 


dT 

~dt 


de 9 , 

Pg r\ H - e G 

dPg 


d(<t>Pg) 


dt 




-dp s 

dt 


- V • (kVT) + V • {<t>p g h 0 y g ) -Q = 0 (33) 


and 


where 


Momentum: {v g ) k + ttk • VP = 0 for k = 1... # dim 

Solid Mass: — rh s = 0 

at 

Gas Mass: + ^ + V ' {<t>P 9 Vg) = 0 

pC v = (1 — P) Pi <C Vv + (3p c C Vc + (j>p g C Vg 


(34) 

(35) 

(36) 

(37) 


and 


and 


e = 


Pv&v PcCc 
Pv Pc 


(38) 


hg (Pg,T) + 



(39) 


Since the model equation for fn s is known, the solid mass conservation equation (Eq. 35) will be substituted 
into the energy and gas mass conservation equations. Consequently, the solution procedure will be to solve for the 
temperature, velocity, and gas density fields through integration of the PDEs given in Eqs. 33, 34, and 36 while the 
solid density will be treated as a nonlinearity, and the resulting ODE governing the solid density evolution will be 
numerically integrated. 


1. Energy Equation 

A Galerkin weak statement can be developed for the energy equation (Eq. 33) by first multiplying it by a suitable 
test function, v, and integrating over the domain O while integrating the 5 th and 6 th terms by parts to give the natural 
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boundary condition terms. 



pC v - P\ 


d<j> de q \ dT 


9 dTdp g ) dt 


+ 


9 

'# dim 


dc 

+ verh s + v ( p g — A + e, 




1 d Pg 


Vv ■ 


d(<t>Pg) 


dt 


dfl 


.k = 1 


(kVr) - Vv ■ (c, bp g h 0 Vg ) 


+ (yh 0g rh w + vq w ) dT = 0 Vv G H ( 


where the boundary (wall) mass flux is 
and the boundary heat flux is 


rn w = (4>Pg) v g ■ n 
q w = ~kVT • n 


vQ 


dil 


(40) 

(41) 

(42) 


2. Momentum Equations 

For the k th momentum equation we again multiply by a suitable test function, v, and integrate over the domain 

f vcj)p ( v g ) k dfl+ I VKk ■ VPdtt = 0 \/v € Hq (43) 

Jn Jo. 

Since the porous flow law is simply a constraint equation, there is no need to integrate by parts to develop natural 
boundary condition terms. While it is recognized that the chain rule could be applied to the pressure gradient term to 
recast it in terms of independent variables, for simplicity it will be assumed that the pressure gradient at the quadrature 
points can be calculated according the local shape function and nodal pressure values. 

3. Gas Mass Conservation Equation 

Likewise, a Galerkin weak statement can be developed for the gas mass conservation equation (Eq. 36). Again, 
the equation will be multiplied by a suitable test function, v, and integrated over the domain f i while integrating the 
3 rd term by parts to give the natural boundary condition term. 

J V ~ ^ V ' ^ + y th w vdr = 0 Vv € Hq (44) 

During the development of CATPISS, it was observed that the hyperbolic nature of the gas continuity equation led to 
numerical instabilities during the solution procedure. Consequently, numerical dissipation or stabilization is required 
to solve the equation system. For the current work, the SUPG (Streamwise Upwind Petrov-Galerkin) stabilization 
technique has been implemented. The primary advantage of using a stabilization technique over simpler numerical 
dissipation techniques is that the solution to the stabilized system is consistent, meaning that the solution to the 
stabilized system still satisfies the original PDEs. Adding the SUPG terms to the gas mass conservation equations 
gives 


( 9{(j>Pg) , 

\ dt ' 


— Vv • ((j)pgV g ) + m s v ) dfl + ® m w vdT 


# elem « 

L / 

e— 1 J n. 


+ E 


TVg • Vv 


(^+m s + V.(# s Vg)) dtt e = 0 Vv € Hq (45) 


where r is the SUPG parameter of intrinsic time scales and is typically based on the element Reynolds number or 
element Peclet number. Expanding the last term gives 

[ ( d{<t>p g ) 

Jn V dt ‘ 


— Vv • ( (f>PgVg ) + rh s v j dfl + ® m w vdY 


# elem « 

c / 

e—1 J fL 




rv g • V v 


( d(<t>Pg) 
\ dt 


+ m s + ((jtpg) V • Vg + Vg • V ( (j> p g ) dCl e = 0 Vv G Hq (46) 
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F. Finite Element Formulation 


Eqs. 40, 43, and 46 can be discretized by expanding the independent variables and test functions in terms of a 
finite dimensional basis 


# nodes 

(x) = (a;) 

7 = 1 

(47) 

# nodes 

0^0)= E V’j (*) 

7 = 1 

(48) 

# nodes 

[K)fc]/ l ( :C ) = E [K)fc]j^(*) 

7 = 1 

(49) 

# nodes 

VP/, (x) = P 3 V ^j( X ) 

7 = 1 

(50) 

# nodes 

V h (x) = ^2 V i^i ( x ) 
i= 1 

(51) 


where the subscript h is introduced to denote a finite dimensional approximation. Since the unknowns are no longer 
functions of x, the PDE system reduces to an ODE system in which the temporal derivatives can be defined as 


# nodes 


r\rji Tr 11UUW 7 

~gf=Yl ( x ) where Tj = — C Tj ) 

3 = 1 


# nodes 


0/7 \ Tr llUU^S 1 

Pa h = E O) where = J t O Pg)j 


dt 


7 = 1 


(52) 

(53) 


d \( V 9)k\h 

dt 


# nodes 

E 

3 = 1 


( V g)k\i^3 ( * ) where [K)j 



(54) 


Since the equation system should be satisfied for all combinations of nodal shape function coefficients, Vi, their choice 
is arbitrary as long as a unique combination is chosen for each node. For the I th nodal equation v, = 1 for i = l and 
Vi = 0 for i jtz l. Consequently, Eqs. 40, 43, and 46 can now become 


# nodes 

E 

3 = 1 


' 2 d<j) de g 

P V Pg dTd P g 

# dim # nodes 


# nodes 


E {<t>Pg)j 


3 = 1 


de 9 , 

P g p b e 0 


ipii/jjdtd 


# nodes 


+ E E [( w fl)fc]j / E T i / 


fc=i j=i 


j'=i 


# nodes 


rr uuucs « C C C 

- E ( < M/)j / h 0g 4>j'Vipi ■ v g dtt + / 'ipiq w dT + / 3j}ih 0g rh w d£ + / ipi [em s - Qj d£l = 0 
^ -j_ « r o 

# nodes „ „ 

E [( v s)fc],- / 4>P' l Pi' t Pjd^+ / 4>iK k ■ VPdfl = 0 

■ =1 J Jo JO 


(55) 


(56) 
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and 


# nodes « # nodes « r r 

E (‘Mt); / V’W’jdO - E / VV’i • ijjjVgdtt + / i/jifn w dT + / 'ip i m s d£l 

Jo Jo Jr Jo 


# nodes # elem 


# nodes # elem 


+ E (^)j E / ^Tv g • + e (M),- E L TVg ■ V0j (iPjV'Vg) dfl e 


# nodes # elem 


+ E E / TVg-Vlpi (Vg • VV’j) dfl e + E / m s TVg-'VipidQ e = 0 


for * = 1,2,...,# nodes. These can more concisely be written ; 


# nodes [* # dim 

E + ('•'d, + E [(«») fc] j M ^ V + T ' K U T + (M)g = 

j=i L /c=i 

# nodes 

E [K)J j K T + F ik = 0 for k = 1 ... # dim 


# nodes 

E [(<MA (E" + E/ sr, ' f; ) + (/C' + K % 

3 - 1 

for * = 1,2,...,# nodes. Where 

^ r cu o. i 


+ if + F p ' SUPG = 


mE = 


/' r 2 ^ 9e g 1 


f t = 


M E = [ p air + e °9 V’iV’jdO 
Jo L u Pg 

M ij V = [ (<f>Pg) ( w s) fe V’iV’jdO 
Jo 

= J VV’i • kV?/yiQ 

= - f KJWi ■ v g dn 

JO 

i pi (ern s - dfl + J tl>iq w dr + J ^h, 

K-f = J fin^i^jdn 

F lk = [ ■ VFdfi 

Jo 

M£ p = [ fafydSl 
Jn 


•, 0a m w dT 


K p - p = - 
*? 


ipj'Vtpi ■ Vgdtt 


A/fP’SU PG = 

ij 


n eiein « 

= E / ^J TV 9 ‘ VV'idOe 


t -p,SUPG _ 


n « # elem « 

/ TVg • V^i OjV-Vg) dO e + E / TVg • V^i (Vg • VV’i) dOg 

♦j r2 e </Q e 
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f p,supg = y I ihsTVs . v ^. dQe 

e—1 J 

F? = / ipiih w dT + / ipirn s dQ, 

Jr J n 


G. Time Discretization 


The semidiscrete weak form of the system given by Eqs. 58 and 60 is discretized in time using backwards fi- 
nite difference schemes. Both first and second-order accurate in time schemes may be derived from Taylor series 
expansions in time about Ft h (t n + 1 ) = Ft n+ 1 : 


FJ n — FJ n+1 


U n - 1 = u r 


dt 

dU n+ 1 
dt 


(tn Wl) + {tn \ n+l)2 + O ((i„ - t n+1 f) 

(tn - i ~ t n+ 1 ) + d2 ^ +1 (t "- 1 - 2 t " +l)2 + O - tn+l) 3 ) 


These expressions can be manipulated as in 2,3 to create difference formulas of the form 


dU n+ 1 
dt 


— OttJJ n + 1 + dtFJ n + ItU n -l + 0 (At£ +1 ) 


to yield either a first or second-order accurate scheme. The weights at, fit, and 7 1 are given for p = 1 and p = 2 in 
Table 1. 

Table 1. First and second-order accurate time discretization coefficients. 

P II a t f3 t 7 1 


^n + l 

- 7t - StTTf + S7 




At n (At n+ i+At„ 


The resulting equation system is 


E ( a t< ; T + K V) T ? +1 + + K Z P ) (# s )" +1 + E [K)*L 


+ E 2? + (#j; + £ AMj* [(»,), 


+ E (#g)" _1 + £ 7tM^ [K) J . + if = 0 = nj (76) 


E [<«.)jr + ^ = ° = for * = 1 - # dim 

i=i 


ff noaes 

E { b ( MP / + Kf UPG ) + + *« SI/PG ] (M); +1 } 

i=i 

# nodes 

+ E {& (j^ij P + MF supg ' s j (4>p g )™ + 7t [fIF p + Mfj SUPG ^ (0p 9 )" -1 } 


+ Ff + F l p ’ SUPG = 0 = (78) 


for * = 1,2,...,# nodes, where 7?., denotes the nonlinear nodal residual equations which is driven to machine zero 
during the iteration process. The integrals in Eqs. 61-74 are evaluated with T n+1 , {(j>p g ) n+l , and v” +1 . 
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H. Linearization 


Through the iteration process the nodal residual equations in Eqs. 76, 77, and 78 will be driven to machine zero so 
that the governing equations are satisfied in a discrete sense. To aid in the iterative process, it is necessary to linearize 
the residual equations in iteration space. This will be done according to the familiar Taylor-series expansion 


K" +1 = W? 


# nodes ( 

E l oRi_ 

I d T, 

3 = 1 l 3 


A Ti 


on, 


0(<t>Pg) 


gjj 


# dim 


A 


on. 


k = 1 




A [(^p) fc ] ■ j + higher order terms 

(79) 


where the superscript v has been introduced to denote iteration level and 


A (), = 0^-0; 


(80) 


Keep in mind that the Jacobian terms, ' , are only nonzero if nodes i and j share an element. 

Since the intent of CATPISS is to allow for easy addition of material models and porous flow laws, the Jacobians 
will not be derived analytically yet they will be calculated exactly (to machine precision) using the complex-step 
method. The implementation of the complex-step method does not disallow the use of analytically derived Jacobian 
terms (or some mix of the two), so these terms can be replaced with analytical expressions at a later date if desired. 

The residuals can be expanded according to a Taylor-series expansion in independent variable space about a point 
[T, ( (f>pg ) , v g ]. For example let’s do the expansion in temperature space and take a complex step , iAT. 


An 1 PPn i B^n 

n[T + iAT, (4>p g ),v g ] = n [T, (# 9 ) , V 9 ]+Z— (AT)— 2 — (AT) 2 - g — (AT) 3 - 


higher order terms 

(81) 


Now looking at just the imaginary part and ignoring the higher order terms gives 

1 7? 

Im {n [T + iAT, (# 9 ) , v a ] } = — . AT - - ^ (AT) 3 (82) 

Solving for the first derivative gives the Jacobian terms 


on 

OT 


Im {n [T + iAT, ( 4>p g ) , v g ]| 
AT 


+ 0 


(A T) z 


(83) 


This process can be repeated for each independent variable. The proper choice for the perturbation steps is not 
immediately obvious from looking at the equations. In general, the perturbation step should have some relationship to 
the order of magnitude of the independent variable. 


A()=r() 


(84) 


Taking advantage of the fact that the problem will be solved on a finite precision machine, r can be chosen so that 
the second order error terms in the Jacobian expressions are smaller the smallest orders represented in the first terms. 
Consequently, the derivatives will be accurate to machine precision. Having such a small step does not affect the divi- 
sion operation in the first term because multiplication/division operations simply result in an exponent shift. Choosing 
r = 10“ 16 for a double precision calculation will accomplish this. 

For each node, the calculation of all of the Jacobian terms will require 

[# Residual Calcs] = [! + (# other nodes on common elements)] x [# d.o.f.] 


complex residual calculations where # dof is the number of degrees of freedom, which is 2 for the current system. 
The scope of the residual recalculations can be reduced by only recalculating those terms that depend on the degree 
of freedom with the current complex perturbation. One fortuitous aspect of this method is that there is no separate 
residual calculation required for the unperturbed state. If only one independent variables is perturbed (use T for 
example), then the real part of Eq. 8 1 can be solved for the residual 


n [T, ( 4>p g ) , v 9 ] = Real {n [T + iAT, ( <j>p g ) , v 9 ]} + O 


(AT) 2 


(85) 
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Again since the independent variable perturbation has been chosen to be sufficiently small, the order of magnitude of 
the error terms is smaller than the smallest order represented in the complex-step residual calculation. This results in 
a residual calculation exact to machine precision. 

The advantage of the complex-step method is in the development cost. The derivation, implementation, and 
debugging of analytical Jacobian terms can be time consuming. Plus, the addition and/or modification of models 
with the complex-step method can be accomplished with relative ease. The drawbacks are the extra time spent doing 
complex arithmetic, and the extra storage required for complex variables. 

I. Boundary Conditions 

CATPISS can handle a rich suite of boundary conditions, including multiple flux type boundary conditions that 
can be summed over a given boundary. CATPISS currently accepts EXODUSII (reference) unstructured grid files 
including side set definitions for boundary conditions. In the current implementation, a side set can be thought of as 
a collection of one or more boundary conditions that is defined for each subdomain exterior boundary. For example, 
a subdomain exterior boundary may have a Dirichlet (specified temperature) condition (1 boundary condition) or 
an Neumann (specified flux) condition ( 1 boundary condition) imposed on it. Alternatively, the boundary could be 
exposed to a convective environment with shock-layer radiation and far-field reradiation (3 boundary conditions) at the 
same time. In addition, CATPISS has the ability to handle both constant and time-dependent boundary conditions each 
with their own input nomenclature. Internally, CATPISS converts constant boundary conditions to time-dependent 
conditions so that their application within the element assembly routines can use the same logic. 


1. Specified Temperature 


Dirichlet conditions and other conditions that require the substitution of a nodal residual equation with a boundary 
condition specific equation will be imposed via the penalty boundary method (PBM) (reference), with this method, 
the L'2 projection of the residuals are added to the matrix. The projection is multiplied by some large factor so that, in 
floating point arithmetic, the existing (smaller) entries in the matrix and right-hand-side are effectively ignored. The 
boundary integral in the weak form becomes 

- [ fi,ndr = o (86) 


e Jr 


where e << 1. For the specified temperature condition, the residual equation is given by 


77 = T spec - T h ( x ) = 0 (87) 

where T spec is the known specified temperature. Substituting in the definition of 7), (a;) in the finite dimensional basis, 
Eq. 47, gives 

# nodes 

77 = T spec - ^2 T j4>j ( x ) = 0 (88) 

i=i 

Substituting the residual definition into Eq. 86 gives the weak form of the discrete residual equations 


77; = 




# nodes 

£*, ( x ) 


dr = o 


(89) 


Employing a numerical integration technique, such as Gaussian quadrature, to evaluate the surface integral gives 


# QPs 


** = 7 £ 


QP= 1 


# nodes 

wfii ( T spec - £ %• 
i= i 


I QP 


(90) 


where QP denotes the quadrature point and to is the integration weight. Linearizing in iteration space according to a 
Taylor-series expansion while ignoring higher order terms gives 


TZ- +1 = 


# nodes 0/T -, 

E dUi 

dP 

3 = 1 3 


A Ti 


( 91 ) 
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Differentiating Eq. 90 gives the Jacobian terms for the linear system 


dUi 

dTj 



( 92 ) 


Given the ease of the Jacobian derivation, this boundary condition has been implemented using analytical Jacobians 
as opposed to using the complex perturbation approach. 


2. Specified Heat Flux 


The specified flux boundary condition is linear and requires no Jacobian contributions. The boundary heat flux 
term in Eq. 66 becomes 

J ipiq w dT = - J ipiq spec dT (93) 


where q spe c is the known specified flux. Note that the sign change is intended to simplify usability. In the original 
PDE derivation, fluxes into the body were considered a negative quantity. So the sign change allows the user to 
input a positive heat flux when the intent is to put heat into the body. Introducing Gaussian quadrature for numerical 
integration gives 


I 'ifriQspecd^ — 


# QPs 

y, -wq apec ipi 
QP—l 


(94) 


3. Convection 


The convective heat flux is given by 


Qconv (T r ec fin ) 


(95) 


Substituting into Eq. 66 gives 


•fiiq w dT 


fiih (T rec 


T w ) dT 


(96) 


Introducing Gaussian quadrature for numerical integration and the finite dimensional basis representation of tempera- 
ture gives 


r # QPs 

/ tpih (T rec - T w ) dF = y 

QP= t 


( # nodes 

Tree - y %•(*) 


I QP 


(97) 


The Jacobian contributions for the convective boundary condition are calculated using the complex perturbation 
method. 


4. ( 

Specified Pressure) Since the independent variables in the linear system are temperature and gas density, the 
specified pressure condition is nonlinear and must be applied through the PBM method to replace the appropriate gas 
density residual equations in the linear system. An additional complication is introduced because the energy equation 
requires a mass flux at the surface to complete the energy balance. Consequently the surface mass flux must be backed 
out from the known data. Before boundary conditions are applied, but after the 

III. Governing Equation Derivation 


A. Energy Equation 

The thermochemical equilibrium energy equation describing the charring porous ablator with pyrolysis gas flow 
on a stationary grid is 

_ V • (kvr) + V • (fip g ho g y g ) -Q = 0 (98) 

where 

pe 0 = (1 — /?) p v e v + Pp c e c + fip g e 0g (99) 
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The virgin and char internal energies are equal to the virgin and char enthalpies respectively and are given by 


D /C = h v/c = h 0 v/c + [ C Vv/c (T')dT' 
Jt° 


and the gas total energy is given by 


e o„ — e g ( pgi T) + 


P is called the “degree of char” or “extent of reaction” and is defined to be 

P = Pv ~ Ps (102) 

Pv Pc 

which represents the char volume fraction within the virgin/char mixture. While there are no distinct virgin and char 
materials in the solid, the extent of reaction is introduced as a modeling assumption to weight thermal properties in 
the decomposition region since only virgin and char properties are known. The extent of reaction definition can be 
rearranged to give the solid density 

p. = (l-p)p v + Pp c (103) 


and the total mixture density is given by 


P — 4>Pg + Ps 


Using the chain rule to evaluate the temporal derivative in the energy equation we have 


d(pe Q ) 


dT , j. de °s , d (^Pg) , - d Ps 


5 T = (i - ® + ^ c - + 


_ PyCy Pc&C 

e = 

Pv Pc 

Concentrating now on the temporal derivative of the gas total energy we have 

de °„ = 

dt dt 9 dt 

where 

deg dT de g dp g 

dt Vg dt d p g dt 

In the solution process, the gas density is a nonlinear function of both T and (<j>p g ) given by 


and consequently 

dp g = dp g d(4>P g ) dp g dT 

dt d (4>p g ) dt dT dt 

Performing the necessary differentiation gives 

9pg = 1 

d(.<t>Pg) 4 1 


dpg = (' <t>pg ) d<t> 

dT 4 > 2 dT 

Substituting into Eq. 107 and collecting terms gives 


de °s _ ( c ^ de a \ dT | 1 de a d (<t>Pg) . 0v g 

dt \ Vg (j) 2 dT dpg ) dt cf) dp g dt 9 dt 
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Substituting back into Eq. 105 gives 


d{pe 0g ) 

dt 


(1 — /?) PvC Vv + (3p c C v c + (4>pg) 



(#g) 9(f) de g Y 
(j> 2 dTdpg) 


O Pg ) de g 
<t> dp g 


dj&Pg) 

dt 


d T 
~dt 


+ {<t>Pg) Vg ■ 


dv g _dp s 


(114) 


Further simplifying gives 


d(pe 


dt 


pC v - p\ 


2 d(j> de g 
9 dfdig\ 


dT 

~dt 


Pg 


deg 

dpg 


dj^Pg) 

dt 


O Pg ) v s 


dv g 

dt 


-dp s 

3 at 


(115) 


where 

pC v = (1 — pi) p u C Vv + (3p c C Vc + 4>pgC Vg (116) 

In summary, the final form of the energy equation is 


PCv - P\ 


2 d4> de g 

dfWgl 


dT 

Ik 


Pg 


deg 

dpg 


d{<t>P g ) 

dt 


O Pg ) v 


dw n 


-dp a 
6 dt 


dt 

- V • (kvr) + V • (<j> Pgho g y g ) - Q = 0 (117) 


where 


and 


pC v — (1 — P) p v C Vv + (3p c C Vc + <j>PgC Vg 

h 0g =hg ( Pg ,T)+^pL 

— Pv&v p c e c 

e = 

Pv Pc 


(118) 

(119) 

( 120 ) 


B. Momentum Equations 

The momentum equations will be simplified through the use of porous flow laws. There are several different 
options, but Darcy’s law will be used for initial development purposes. 

1. Darcy’s Law 

Darcy’s law was initially empirically developed through porous flow experiments by Darcy, and it has since been 
analytically derived directly from the Navier-Stokes equations. Darcy’s law governs the volumetric flow rate of a 
steady incompressible low Reynolds number laminar fluid flow through a uniform porous medium. While it is rec- 
ognized that the pyrolysis gas flow through the char layer is neither steady nor incompressible, Darcy’s law is the 
simplest porous flow law and will be used for initial development. More appropriate momentum equation models can 
be implemented in the future. 

The volumetric flow rate, Q, is given by 

Q = -A — VP (121) 

P 

Consequently, the mass flow rate is 

PgQ = ~ Pg A-VP (122) 

p 

Averaging the mass flow rate over the sample are through which is flows gives the mass flux 

pg<, = -Pg K ^ P (123) 

P 
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where \' g is the superficial or filtration velocity of the gas. If it is assumed that the surface porosity is equal to the 
volumetric porosity, then the average velocity (or pore velocity) of the gas is given by 



The mass flux can now be restated in terms of the gas velocity 


(124) 


= ~pa-y p 

h 1 


There will be one equation in the PDE system for each velocity component 

( v g)k + *k ' VP = 0 

where k represents the Cartesian dimension and is the k th row vector of the permeability tensor 


k(/3) = 


Kn (/3) «i2 (/?) k 13 (/3) 

«21 (/3) K 22 (/?) H23 (/3) 

«3t (/3) «32 (/?) K 33 (/3) 


where each component is a function of the extent of reaction. 

C. Gas Mass Conservation Equation 


d<t>Pg 

dt 


= -rn s 


V • (#gVg) 


IV. Verification 


(125) 


(126) 


(127) 


(128) 


Multiple levels of verification are performed to ensure that CATPISS correctly solves the governing equations 
over a variety of geometries and flow conditions. First, temperature profiles are obtained for analytic solutions to 
the one-dimensional heat equation for several different boundary conditions. Next, multidimensional solutions are 
calculated for the convective boundary condition on a cube, cylinder and sphere then compared with the temperature 
distributions estimated with CATPISS. Second-order convergence of the spatial discretization routine in CATPISS 
is demonstrated in the normalized error between the exact solution and code output as the grid is refined. Lastly, 
the method of manufactured solutions provides general evidence that the code accurately recreates a preconceived 
temperature profile as defined by a source term in the governing equation. 


A. One-dimensional Heat Equation 

The one -dimensional heat equation describes the distribution and evolution of temperature in a medium 


ar 

at ax 1 


0 < x < L 


(129) 


where p is the density of the material, C v is the specific heat at constant volume and k is the thermal conductivity. 
The heat equation is of particular interest in this study since it is in the same form as the ablation governing equa- 
tions and has exact solutions for a myriad of boundary conditions. A comparison of the temperature profile given by 
the exact solution and simulated numerically by CATPISS is presented for a quasi two-dimensional surface assum- 
ing unsteady conduction with constant thermal material properties, unsteady conduction with temperature-dependent 
material properties, transient conduction with standard wall convection, transient conduction with a specified wall 
temperature, steady-state radiation and reradiation. Exact solutions also exist for simple three-dimensional geometries 
with exclusively adiabatic and convective surface boundary conditions. Although the analytic solutions are provided 
in this report and discussed briefly for each case, the derivation and implementation details are omitted for brevity 
except where numerical error may be a notable factor. 
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B. Unsteady Conduction with Constant Properties 


The unsteady conduction problem requires Robin boundary conditions to specify the heat flux each side of the 
two-dimensional plate. The side defined as x = 0 is arbitrarily selected as the surface to which there is a positive heat 
flux while the remaining three sides are assumed to be adiabatic as described by the following expressions: 


dT 

i i •// 

- k dx\*=°- q 

(130) 

dT 

~ k ~dx \y=°’ x=L 'y= L = ^ 

The initial condition can be set to any uniform temperature such that 

(131) 

E-? 

II 

o' 

II 

h 

(132) 


The exact solution in non-dimensional form is 


T (x, t) — Tq at 1 x 
q" L/k ~ L 2 + 3 ~ £ 



(133) 


where T 0 is the initial temperature, L is the characteristic length of the plate and a represents the thermal diffusivity of 
the material. The general procedure for calculating an infinite sum in any of the exact solutions is to continue adding 
terms until the absolute difference between consecutive values is within machine precision or IE-15. The material 
properties and other input parameters for the unsteady conduction problem are provided in Table ?? and are consistent 
for each of the cases presented in this report. 


C. Unsteady Conduction with Variable Properties 

While the previous example was derived based on constant material properties, it is also possible to determine an 
exact solution wherein the specific heat and thermal conductivity are allowed to vary as a function of temperature to 
render a more physically realistic result. Despite evidence suggesting that both specific heat and thermal conductivity 
may actually have a quadratic thermal dependence, it is simply assumed for the purposes of this investigation that the 
material properties exhibit a linear relationship to temperature. Either way, CATPISS permits the user to specify each 
of the input material properties at any number of temperatures to achieve the desired level of accuracy. 

k(T) = k 1 + ^^(T-T 1 ) 

C v C T ) = C v , t + ° V T 2 ~ i?’ 1 (T - Ti) 

_ ki _ fc 2 _ k (T) 
pCy^l pC Vt 2 pC v (T) 

The one -dimensional heat equation for temperature-dependent material properties is 

T d ( dT\ 

pC ^ = dx\ k dx) °^ L (13?) 


(134) 

(135) 

(136) 


First, Kirchoff's transformation is applied to define a new variable 9: 

9 = - — f k (f) df (138) 

^ ref J T re f ^ 

In a similar manner to the previous case, both heat flux boundary conditions and an initial temperature profile must be 
specified prior to searching for an exact solution: 

-A ;^\ x=0 = g" (139) 

dO 

— k — \ y =o,x=L,y=L = 0 (140) 
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( 141 ) 


0( x ,t = 0) = d o = 0(T = T o ) 

Several intermediate steps are taken to manipulate the original heat equation eventually leading to the result 


6 (x, t) - 6> 0 
q"L/k 


at 

Z 2 


x 

L 


& 


1 


■ exp 


2 2 
—nn 


( x \ 

cos ^n7r — ) 


(142) 


It is interesting to note that this solution is in fact the same as for the constant property unsteady conduction problem 
with T replaced by 9. Naturally, however, the results are quite different and deserve special consideration. As time 
advances the two cases tend to diverge since the medium’s temperature is continually increasing. 

D. Transient Conduction with Standard Convection 

The problem of transient conduction with standard convection involves Robin boundary conditions such that one 
side of the plate is exposed to a fluid of temperature T r and the remaining three sides are adiabatic: 


-fc^|*=o = H[T(t)-T r ] 
dT 

~ k~Q^\ y=0,x=L,y=L = 0 


(143) 


(144) 


Once again, the initial condition is uniform throughout the plate and set to a value denoted as 7 1,. The exact solution is 


T ( x , t) — T r 
Tn - T r 




n= 1 


sm v n 


v n + sin v n cos v n 


- vied 

exp | — — | cos 


L 2 


m 


where v n is the index eigenvalue that satisfies the equation 

v„ sin Vn = Bi cos Vn 


(145) 


(146) 


The parameter Bi in Eq. 146 is known as the Biot number and is defined as the ratio hL/k. In order to evaluate the 
terms of the infinite sum, it is first necessary to find the value of v n that completes the implicit relationship for each 
index n. In general, this type of parameter is calculated using the Newton-Raphson method with an initial guess of 
Bi-7r/2 and approximate periodicity of 7 r. Implementation of the numerical technique does not appear to noticeably 
affect the efficiency of the exact solution routine which is written in the FORTRAN programming language. Moreover, 
the convergence tolerance between successive root estimations is set sufficiently high to avoid the propagation of error 
in the infinite sum. 

E. Transient Conduction with Specified Temperature 

The exact solution for transient conduction with a specified wall temperature is similar in form to the standard 
convection case although cast in terms of a freestream instead of recovery temperature 7’,n/ and somewhat simpler. 
The boundary and initial conditions are set as 


resulting in an analytic solution 

T (x, t) - T„ 


dT 

^ dx ^ x — T °° 
dT 

~ k-Qy\y=°’ x = L ’y= L ~ 0 
T (x, t = 0) = ^ 

(-l ) n+1 


To — 77, 


= 2£ 


n—1 


6XP ) C ° S l 


v n x 

~L~ 


where v n is now the index eigenvalue of the equation 


cos v n = 0 


v n = (2 n - 1) - 


(147) 

(148) 

(149) 

(150) 

(151) 

(152) 


The computation cost is reduced since the value of v n is automatically known for all possible indices n without the 
need of numerical approximation. 
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F. Steady-State Radiation and Reradiation 


Both steady-state radiation and reradiation are governed by the total differential equation 


dPT 

dx 2 


= 0 


For the case of radiation, the constant heat flux q is directed into the medium with an absorptivity a: 


(153) 


Qrad = aq (154) 

Reradiation implies that constant heat flux is exiting the medium with an emissivity e and background sink temperature 
1), in addition to two parallel adiabatic sides: 

- fc S l "= 0 = £CT (' T res ~ T f ) 

.dr, _ T 
k dJ X=L ~ 6 
_,dT _ n 

k dx ' y= ®' y=L 0 

The reradiaton exact solution, therefore, is simply 

ea (T 4 s - Tf) = j{T f - T b ) (158) 

which is found via the Newton-Raphson approach to machine tolerance. A so-called unsteady reradiation verification 
case can be forged by extracting the evolution of surface temperature for unsteady conduction and assuming a linear 
relationship between emissivity and temperature to calculate the time history of sink temperature. 


(155) 

(156) 

(157) 


G. Three-Dimensional Heat Equation on a Cube 

In the case of transient conduction with standard convection, it is possible to extend the one-dimensional solution 
to a cubic geometry through the principle of superposition. The exact solution given in Eq. 145 is implemented in the 
x, y and z directions independently to render the functions f(x), f(y) and f(z) which are then multiplied to give the 
non-dimensional solution at any point in the cube: 

T{x y*i- Tr =/(*>•/(,)■/(») cm 

It is also possible to specify different properties and heat flux conditions on each face of the cube, but for the purposes 
of this study opposite faces will be identical such that the center plane acts as an adiabatic surface to mimic the back 
of an ablating thermal protection material. The characteristic length L and convection coefficient h are permitted to 
vary between each direction. 


H. Three-Dimensional Heat Equation on a Cylinder 


Transient conduction with standard convection on a cylinder is considerably more complicated in both the nature 
and implementation of the exact solution. First, several dimensionless parameters are defined 


r 

V 


hr 0 at 

Bl =^T’ Fo =Z2 


(160) 


where r is the cylinder radius, Bi is the Biot number as described previously and Fo is the Fourier number. Transform- 
ing the one -dimensional heat equation into cylindrical coordinates renders a Bessel-type partial differential equation 
with the solution 


T (r, t) - Tqq 
Ti - ^ 


OO 

E 

n= 


] K n Jo ( K bn ~^j eXP (~ b n Fo ) 


(161) 


As is the case with convection in one-dimension, b n and K n are eigenvalue index parameters but now depend on 
Bessel functions of the first kind J: 


bnJi (pn) — BiJ 0 (b n ) 


(162) 
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( 163 ) 


" (bl + BP) J 0 (b n ) 

It should be noted that this analysis is only valid for long cylinders, or cylinders in which the length is much larger 
than the radius. Finite cylinders require an intersection of Eq. 161 with the non-dimensional temperature profile for a 
convective plate of half the characteristic length L (Eq. 159) as performed in the following manner: 

6 (r,x,t) = T (r,x,t) - Too (164) 

9 t = T i -T oa (165) 


'9 (r,x, t)~ 

= 

’ 9(r , t)' 


9 (x, t) 

(166) 


finite cyl (L) 


long cyl 

0i 

plate (L/2) 


The main difficulty with this approach is calculating the value of various order Bessel functions at every point 
within the domain in addition to determining the eigenvalue index values with an inherently iterative numerical tech- 
nique. One option commonly employed is to tabulate regular intervals of the Bessel functions and interpolate to the 
desired location. Although much more efficient than converging the infinite series definition of Bessel functions of the 
first kind, such a method greatly limits the applicability of the verification software to geometries of a similar size and 
simplicity. An optimized parallel processor code is a particularly attractive resolution to this issue. 
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